function [smooth_st,etamat,stt,countmat,handvec,flag_mode,topstr,countobs]=...
    mode_countercont(funcmod,fullp,Y,stloc,stanames,shocell,sample,outfold,flag_mode,flag_noplot)
% =========================================================================
%function [smooth_st,etamat,stt,countmat,handvec,flag_mode,topstr,countobs]=...
%    mode_countercont(funcmod,fullp,Y,stloc,stanames,shocell,sample,outfold,flag_mode,flag_noplot)
%
% FLAG_MODE determines how to do the counterfactuals. Mode 1 is default 
% If OUTFOLD not passed, graphs will not be saved 
%
% 1.A Solve model at fullp and obtain states and innovation using Kalman Filter 
% Then compute counterfactual paths of observables and states as either 
% Shut down one shock a time   (FLAG_MODE=1, default) 
% Shut all      shocks but one (FLAG_MODE=2         )
% COUNTMAT      [nobs nst nx]   reports the counterfactual for those states
%               selected with row indicator STLOC. 
%               added to the code output. 
% HADNVEC       Vector of figure handles 
% FLAG_MODE     If 1 or 2 
% TOPSTR        Sting Indicating Mode 
% COUNTOBS is the counterfactual for the observables
%
% October 23 2007. Modified KFILTER_FULL.m to ensure counterfactual
% observables and states add up to actual data & smoothed states 
% Note: This will be true only if FLAG_MODE == 1 
% 
% Warning: Counterfactual on Y is on demeaned data to check that these are recovered 
% correctly (when summed across shocks). Hence if wish to form levels from counterfactuals 
% and the model has constants in C, the resulting levels will correspond to linearly 
% DETRENDED observables. 
%
% Alejandro Justiniano October 23 2007 
% March 2010: produce flag_nond and warning if SDX is non-diagonal; will
% not recover exact decomposition of states. ADDED FLAG_NOPLOT
% =======================================================================
cucd=pwd; 
if nargin < 10 || isempty(flag_noplot) 
    flag_noplot=0; 
end
if nargin < 9 || isempty(flag_mode);
    flag_mode=1;
elseif flag_mode~=1
    flag_mode=2;
end
if flag_mode==2;
    disp('Exclude all but one shock');pause(0.2);topstr=[' is excluded']; 
    quer('c'); 
else 
    disp('One shock on at a time');pause(0.2);topstr=[' only is on']; 
end 
topstr=[' counterfactual, shock(i) ',topstr,'    actual (dark), counter (gray)']; 
if nargin < 8 || isempty(outfold); 
    flag_save=0; 
else 
    cd(outfold);cd(cucd); flag_save=1; 
end 
[GG,CC,RR,eu,SDX,ZZ]=play_eugen(funcmod,fullp);
if ~isequal(eu,[1;1]);error('Model does not have unique stable solution');end
nst=length(stloc);nobs=size(Y,1);nx=size(SDX,1); 
if length(shocell)~=nx;error('Mistmatch SHOCELL and NX');end
if length(stanames)~=nst;error('Mistmatch stanames and STLOC');end
ydem=Y-repmat( (ZZ*CC)',[nobs 1] );
[stt,etamat,smooth_st,junk,junk,junk,countszero]=kfilter_full(ydem,ZZ,GG,RR,SDX,[]);
clear junk; 
countmat=zeros(nobs,nst,nx); 
countobs  =zeros(nobs,size(ZZ,1),nx); 
% Begin Counterfactual 
for ii=1:nx 
    tempeta=zeros(nx,nobs);
    if flag_mode==2
        tempind=(1:nx);tempind(ii)=[];
        tempeta(tempind,:)=etamat(:,tempind)';
    else
        tempeta(ii,:)=etamat(:,ii)';
    end
    % ==================================================
    % Loop follows KFILTER_FULL.m 
    ysim=zeros(size(ZZ,1),nobs);
    stsim=zeros(size(GG,1),nobs); 
    stsim(:,1)=countszero(:,ii); 
    ysim(:,1)=ZZ*stsim(:,1);
    for jj=2:nobs;
        stsim(:,jj)=GG*stsim(:,jj-1) + RR*(tempeta(:,jj));
        ysim(:,jj)=ZZ*stsim(:,jj);
    end
    % =====================================================
    countmat(:,:,ii)=stsim(stloc,:)'; 
    countobs(:,:,ii)=ysim'; 
end 
% Check that the sum of the counterfactual matrices adds up to the
% smoothed states October 23 2007 
% If FLAG_MODE == 2 the counterfactual need not recover the observables 

% Check if there are non-diagonal elements in SDX. In this case will NOT
% recover an exact decomposition of shocks 
if ~isempty( find( triu(SDX,1)~=0 ) )
    disp('Warning SDX has non-diagonal elements, will not recover exact decomposition of states'); 
    pause(1); 
    flag_nond=1; 
else 
    flag_nond=0; 
end


if flag_mode==1
    checkysim=max(abs(sum(countobs,3)-ydem));
    dispaj('Max difference counterfactual & observables is:',max(checkysim));
    if max(checkysim) > 1e-4;
        if flag_nond==0
            error('Cannot recover original states with counterfactual');
        else
            disp('Cannot recover original states with counterfactual as errors are correlated');
        end
        
    end
    checksim=max(abs(sum(countmat,3)-smooth_st(:,stloc)));
    dispaj('Max difference counterfactual & smoothed states:',max(checksim));
    if max(checkysim) > 1e-4;
        if flag_nond==0
            error('Cannot recover original states with counterfactual');
        else
            disp('Cannot recover original states with counterfactual as errors are correlated');
        end
        
    end

else
    disp('FLAG_MODE==2 hence will not recover observables')
end
% Plot actual and counterfactual 
if flag_noplot==1; 
    disp('FLAG_NOPLOT==1, no plots will be produced'); 
    handvec=[]; 
    return 
end 

ncol=2; 
nrow=ceil(nx/2); 
handvec=zeros(nst,1); 
for ii=1:nst 
    figure;
    for jj=1:nx   
        subplot(nrow,ncol,jj);
        hand=plot(sample,smooth_st(:,stloc(ii)),'LineWidth',2);hold on; 
        set(hand,'Color',[0 0 0])
        hand=plot(sample,squeeze(countmat(:,ii,jj)),'LineWidth',2.5); 
        set(hand,'Color',[0.7 0.7 0.7]);  hold off; 
        %set(gca,'XGrid','on'); 
        title(shocell{jj}); 
        axis tight; 
    end 
    hand=suptitle_withpatch([stanames{ii},topstr]);
    set(hand,'FontSize',12); 
    handvec(ii)=gcf;
    set(handvec(ii),'PaperPosition',[0.75 0.5 7.5 10])
    if flag_save==1
        cd(outfold);
        print(gcf,'-dpdf',['counter ',stanames{ii}]);
        cd(cucd);
    end
end 